Data and code for paper on evolution during a 2017 epidemic of Pasteuria ramosa in Little Appleton Lake. Authors: Camden D. Gowler Haley Essington Patrick A. Clay Bruce O’Brien Clara L. Shaw Rebecca Bilich Meghan A. Duffy
lilApp <- readr::read_csv(here("data/lilApp.csv"))
## New names:
## * `` -> ...1
## Rows: 11 Columns: 54
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (3): Date, Lake, Counted.By
## dbl (51): ...1, Total, X, UA, UJ, Uninfected.Males, Uninfected.Ephip, A.Past...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# do this to see which parasites are there (past, micg, spiro)
ttt <- lilApp %>%
select(Uninfected.Males:M.Past...Spider) %>%
gather("Parasite", "Count",Uninfected.Males:M.Past...Spider) %>%
group_by(Parasite) %>%
summarise(totinfe = sum(Count)) %>%
filter(totinfe > 0)
past <- lilApp %>% # find columns with pasteuria in title
select(contains("past"))
lilApp$pasteuria.inf <- rowSums(past)
micg <- lilApp %>%
select(contains("mic"))
lilApp$micg.inf <- rowSums(micg)
spiro <- lilApp %>%
select(contains("spiro"))
lilApp$spiro.inf <- rowSums(spiro)
LilApp_allpara <- lilApp %>%
mutate(pasteuria.prev = pasteuria.inf/Total,
micg.prev = micg.inf/Total,
spiro.prev = spiro.inf/Total)
# super contrived way of doing this
sample_vec <- c(rep(NA, 11))
collect <- c(3, 5, 7)
for (i in 1:length(collect)) {
x <- collect[i]
sample_vec[x] <- paste0("collected",i)
}
sample_vec[is.na(sample_vec)] <- "regular"
LilApp_allpara$sampling <- sample_vec
LilApp_allpara$sampling <- as.factor(LilApp_allpara$sampling)
str(LilApp_allpara)
## spec_tbl_df [11 × 61] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ ...1 : num [1:11] 1 2 3 4 5 6 7 8 9 10 ...
## $ Date : chr [1:11] "7/27/17" "8/10/17" "8/28/17" "9/8/17" ...
## $ Lake : chr [1:11] "LittleAppleton" "LittleAppleton" "LittleAppleton" "LittleAppleton" ...
## $ Counted.By : chr [1:11] "MAR" "KKH" "MAR" "KKH" ...
## $ Total : num [1:11] 265 308 359 321 206 249 116 171 203 154 ...
## $ X : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
## $ UA : num [1:11] 70 97 49 46 22 101 36 62 90 41 ...
## $ UJ : num [1:11] 185 179 183 163 71 51 30 83 100 88 ...
## $ Uninfected.Males : num [1:11] 0 0 0 0 0 2 1 13 12 22 ...
## $ Uninfected.Ephip : num [1:11] 0 0 0 0 0 0 0 0 0 3 ...
## $ A.Pasteuria : num [1:11] 0 7 6 17 18 3 5 3 1 0 ...
## $ J.Pasteuria : num [1:11] 0 10 61 82 55 79 31 5 0 0 ...
## $ Male.Pasteuria : num [1:11] 0 0 0 0 0 1 0 0 0 0 ...
## $ Ephip.Pastueria : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
## $ Brood : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
## $ E.Brood : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
## $ A.MicG : num [1:11] 8 13 39 12 19 11 10 5 0 0 ...
## $ J.MicG : num [1:11] 0 0 20 1 13 1 0 0 0 0 ...
## $ Male.MicG : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
## $ Ephip.MicG : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
## $ A.Spiro : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
## $ J.Spiro : num [1:11] 2 2 1 0 0 0 0 0 0 0 ...
## $ Male.Spiro : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
## $ Ephip.Spiro : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
## $ A.White.Spiro : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
## $ J.White.Spiro : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
## $ Male.White.Spiro : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
## $ E.White.Spiro : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
## $ A.Yellow.Spiro : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
## $ A.Metsch : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
## $ J.Metsch : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
## $ Male.Metsch : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
## $ Ephip.Metsch : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
## $ A.Larssonia : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
## $ A.Spiderman : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
## $ J.Spiderman : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
## $ E.Spiderman : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
## $ Male.Spiderman : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
## $ A.MicG...Past : num [1:11] 0 0 0 0 4 0 1 0 0 0 ...
## $ J.MicG...Past : num [1:11] 0 0 0 0 4 0 2 0 0 0 ...
## $ A.MicG...Spiro : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
## $ A.Metsch...Brood : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
## $ AMicG...Brood : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
## $ A.MicG...Metsch : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
## $ J.MicG...Metsch : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
## $ M.MicG...Metsch : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
## $ A.Metsch...Past : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
## $ J.Metsch...Past : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
## $ A.MicG...White.Spiro : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
## $ Male.white.spiro...metsch: num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
## $ A.Past...Spider : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
## $ M.Past...Spider : num [1:11] 0 0 0 0 0 0 0 0 0 0 ...
## $ Julian.Day : num [1:11] 208 222 240 251 261 276 280 286 300 304 ...
## $ Density : num [1:11] 235212 305732 107071 128571 67081 ...
## $ pasteuria.inf : num [1:11] 0 17 67 99 81 83 39 8 1 0 ...
## $ micg.inf : num [1:11] 8 13 59 13 40 12 13 5 0 0 ...
## $ spiro.inf : num [1:11] 2 2 1 0 0 0 0 0 0 0 ...
## $ pasteuria.prev : num [1:11] 0 0.0552 0.1866 0.3084 0.3932 ...
## $ micg.prev : num [1:11] 0.0302 0.0422 0.1643 0.0405 0.1942 ...
## $ spiro.prev : num [1:11] 0.00755 0.00649 0.00279 0 0 ...
## $ sampling : Factor w/ 4 levels "collected1","collected2",..: 4 4 1 4 2 4 3 4 4 4 ...
## - attr(*, "spec")=
## .. cols(
## .. ...1 = col_double(),
## .. Date = col_character(),
## .. Lake = col_character(),
## .. Counted.By = col_character(),
## .. Total = col_double(),
## .. X = col_double(),
## .. UA = col_double(),
## .. UJ = col_double(),
## .. Uninfected.Males = col_double(),
## .. Uninfected.Ephip = col_double(),
## .. A.Pasteuria = col_double(),
## .. J.Pasteuria = col_double(),
## .. Male.Pasteuria = col_double(),
## .. Ephip.Pastueria = col_double(),
## .. Brood = col_double(),
## .. E.Brood = col_double(),
## .. A.MicG = col_double(),
## .. J.MicG = col_double(),
## .. Male.MicG = col_double(),
## .. Ephip.MicG = col_double(),
## .. A.Spiro = col_double(),
## .. J.Spiro = col_double(),
## .. Male.Spiro = col_double(),
## .. Ephip.Spiro = col_double(),
## .. A.White.Spiro = col_double(),
## .. J.White.Spiro = col_double(),
## .. Male.White.Spiro = col_double(),
## .. E.White.Spiro = col_double(),
## .. A.Yellow.Spiro = col_double(),
## .. A.Metsch = col_double(),
## .. J.Metsch = col_double(),
## .. Male.Metsch = col_double(),
## .. Ephip.Metsch = col_double(),
## .. A.Larssonia = col_double(),
## .. A.Spiderman = col_double(),
## .. J.Spiderman = col_double(),
## .. E.Spiderman = col_double(),
## .. Male.Spiderman = col_double(),
## .. A.MicG...Past = col_double(),
## .. J.MicG...Past = col_double(),
## .. A.MicG...Spiro = col_double(),
## .. A.Metsch...Brood = col_double(),
## .. AMicG...Brood = col_double(),
## .. A.MicG...Metsch = col_double(),
## .. J.MicG...Metsch = col_double(),
## .. M.MicG...Metsch = col_double(),
## .. A.Metsch...Past = col_double(),
## .. J.Metsch...Past = col_double(),
## .. A.MicG...White.Spiro = col_double(),
## .. Male.white.spiro...metsch = col_double(),
## .. A.Past...Spider = col_double(),
## .. M.Past...Spider = col_double(),
## .. Julian.Day = col_double(),
## .. Density = col_double()
## .. )
## - attr(*, "problems")=<externalptr>
summary(LilApp_allpara)
## ...1 Date Lake Counted.By
## Min. : 1.0 Length:11 Length:11 Length:11
## 1st Qu.: 3.5 Class :character Class :character Class :character
## Median : 6.0 Mode :character Mode :character Mode :character
## Mean : 6.0
## 3rd Qu.: 8.5
## Max. :11.0
## Total X UA UJ Uninfected.Males
## Min. :116.0 Min. :0 Min. : 22.00 Min. : 30.0 Min. : 0.000
## 1st Qu.:187.0 1st Qu.:0 1st Qu.: 43.50 1st Qu.: 75.5 1st Qu.: 0.000
## Median :220.0 Median :0 Median : 62.00 Median : 88.0 Median : 1.000
## Mean :233.8 Mean :0 Mean : 62.73 Mean :110.3 Mean : 6.636
## 3rd Qu.:286.5 3rd Qu.:0 3rd Qu.: 83.00 3rd Qu.:171.0 3rd Qu.:12.500
## Max. :359.0 Max. :0 Max. :101.00 Max. :185.0 Max. :23.000
## Uninfected.Ephip A.Pasteuria J.Pasteuria Male.Pasteuria
## Min. : 0.000 Min. : 0.000 Min. : 0.00 Min. :0.00000
## 1st Qu.: 0.000 1st Qu.: 0.500 1st Qu.: 0.00 1st Qu.:0.00000
## Median : 0.000 Median : 3.000 Median :10.00 Median :0.00000
## Mean : 3.909 Mean : 5.455 Mean :29.36 Mean :0.09091
## 3rd Qu.: 0.000 3rd Qu.: 6.500 3rd Qu.:58.00 3rd Qu.:0.00000
## Max. :40.000 Max. :18.000 Max. :82.00 Max. :1.00000
## Ephip.Pastueria Brood E.Brood A.MicG J.MicG
## Min. :0 Min. :0 Min. :0 Min. : 0.00 Min. : 0.000
## 1st Qu.:0 1st Qu.:0 1st Qu.:0 1st Qu.: 3.00 1st Qu.: 0.000
## Median :0 Median :0 Median :0 Median :10.00 Median : 0.000
## Mean :0 Mean :0 Mean :0 Mean :10.73 Mean : 3.182
## 3rd Qu.:0 3rd Qu.:0 3rd Qu.:0 3rd Qu.:12.50 3rd Qu.: 1.000
## Max. :0 Max. :0 Max. :0 Max. :39.00 Max. :20.000
## Male.MicG Ephip.MicG A.Spiro J.Spiro Male.Spiro
## Min. :0 Min. :0 Min. :0 Min. :0.0000 Min. :0
## 1st Qu.:0 1st Qu.:0 1st Qu.:0 1st Qu.:0.0000 1st Qu.:0
## Median :0 Median :0 Median :0 Median :0.0000 Median :0
## Mean :0 Mean :0 Mean :0 Mean :0.4545 Mean :0
## 3rd Qu.:0 3rd Qu.:0 3rd Qu.:0 3rd Qu.:0.5000 3rd Qu.:0
## Max. :0 Max. :0 Max. :0 Max. :2.0000 Max. :0
## Ephip.Spiro A.White.Spiro J.White.Spiro Male.White.Spiro E.White.Spiro
## Min. :0 Min. :0 Min. :0 Min. :0 Min. :0
## 1st Qu.:0 1st Qu.:0 1st Qu.:0 1st Qu.:0 1st Qu.:0
## Median :0 Median :0 Median :0 Median :0 Median :0
## Mean :0 Mean :0 Mean :0 Mean :0 Mean :0
## 3rd Qu.:0 3rd Qu.:0 3rd Qu.:0 3rd Qu.:0 3rd Qu.:0
## Max. :0 Max. :0 Max. :0 Max. :0 Max. :0
## A.Yellow.Spiro A.Metsch J.Metsch Male.Metsch Ephip.Metsch A.Larssonia
## Min. :0 Min. :0 Min. :0 Min. :0 Min. :0 Min. :0
## 1st Qu.:0 1st Qu.:0 1st Qu.:0 1st Qu.:0 1st Qu.:0 1st Qu.:0
## Median :0 Median :0 Median :0 Median :0 Median :0 Median :0
## Mean :0 Mean :0 Mean :0 Mean :0 Mean :0 Mean :0
## 3rd Qu.:0 3rd Qu.:0 3rd Qu.:0 3rd Qu.:0 3rd Qu.:0 3rd Qu.:0
## Max. :0 Max. :0 Max. :0 Max. :0 Max. :0 Max. :0
## A.Spiderman J.Spiderman E.Spiderman Male.Spiderman A.MicG...Past
## Min. :0 Min. :0 Min. :0 Min. :0 Min. :0.0000
## 1st Qu.:0 1st Qu.:0 1st Qu.:0 1st Qu.:0 1st Qu.:0.0000
## Median :0 Median :0 Median :0 Median :0 Median :0.0000
## Mean :0 Mean :0 Mean :0 Mean :0 Mean :0.4545
## 3rd Qu.:0 3rd Qu.:0 3rd Qu.:0 3rd Qu.:0 3rd Qu.:0.0000
## Max. :0 Max. :0 Max. :0 Max. :0 Max. :4.0000
## J.MicG...Past A.MicG...Spiro A.Metsch...Brood AMicG...Brood A.MicG...Metsch
## Min. :0.0000 Min. :0 Min. :0 Min. :0 Min. :0
## 1st Qu.:0.0000 1st Qu.:0 1st Qu.:0 1st Qu.:0 1st Qu.:0
## Median :0.0000 Median :0 Median :0 Median :0 Median :0
## Mean :0.5455 Mean :0 Mean :0 Mean :0 Mean :0
## 3rd Qu.:0.0000 3rd Qu.:0 3rd Qu.:0 3rd Qu.:0 3rd Qu.:0
## Max. :4.0000 Max. :0 Max. :0 Max. :0 Max. :0
## J.MicG...Metsch M.MicG...Metsch A.Metsch...Past J.Metsch...Past
## Min. :0 Min. :0 Min. :0 Min. :0
## 1st Qu.:0 1st Qu.:0 1st Qu.:0 1st Qu.:0
## Median :0 Median :0 Median :0 Median :0
## Mean :0 Mean :0 Mean :0 Mean :0
## 3rd Qu.:0 3rd Qu.:0 3rd Qu.:0 3rd Qu.:0
## Max. :0 Max. :0 Max. :0 Max. :0
## A.MicG...White.Spiro Male.white.spiro...metsch A.Past...Spider M.Past...Spider
## Min. :0 Min. :0 Min. :0 Min. :0
## 1st Qu.:0 1st Qu.:0 1st Qu.:0 1st Qu.:0
## Median :0 Median :0 Median :0 Median :0
## Mean :0 Mean :0 Mean :0 Mean :0
## 3rd Qu.:0 3rd Qu.:0 3rd Qu.:0 3rd Qu.:0
## Max. :0 Max. :0 Max. :0 Max. :0
## Julian.Day Density pasteuria.inf micg.inf
## Min. :208.0 Min. : 4013 Min. : 0.00 Min. : 0.00
## 1st Qu.:245.5 1st Qu.: 8457 1st Qu.: 0.50 1st Qu.: 3.00
## Median :276.0 Median : 14333 Median :17.00 Median :12.00
## Mean :267.7 Mean : 81362 Mean :35.91 Mean :14.91
## 3rd Qu.:293.0 3rd Qu.:117821 3rd Qu.:74.00 3rd Qu.:13.00
## Max. :317.0 Max. :305732 Max. :99.00 Max. :59.00
## spiro.inf pasteuria.prev micg.prev spiro.prev
## Min. :0.0000 Min. :0.000000 Min. :0.00000 Min. :0.000000
## 1st Qu.:0.0000 1st Qu.:0.002463 1st Qu.:0.01689 1st Qu.:0.000000
## Median :0.0000 Median :0.055195 Median :0.04050 Median :0.000000
## Mean :0.4545 Mean :0.151335 Mean :0.06050 Mean :0.001530
## 3rd Qu.:0.5000 3rd Qu.:0.320872 3rd Qu.:0.08013 3rd Qu.:0.001393
## Max. :2.0000 Max. :0.393204 Max. :0.19417 Max. :0.007547
## sampling
## collected1:1
## collected2:1
## collected3:1
## regular :8
##
##
LA_pastprev_plot <- ggplot() +
geom_line(data = LilApp_allpara, aes(x = Julian.Day, y = pasteuria.prev*100), size = 1.1, alpha = 0.7) +
geom_point(data = LilApp_allpara, aes(x = Julian.Day, y = pasteuria.prev*100, fill=sampling), size = 5, alpha = 1.0, shape = 21, show.legend = FALSE) +
ylim(0, 0.45*100) +
scale_x_continuous(breaks = c(213, 244, 274, 305), labels = c("Aug-1", "Sept-1","Oct-1", "Nov-1")) +
scale_fill_manual(values=c("#5445b1", "#749dae", "#f3c483",'#7E7E7E')) +
xlab("") +
ylab("Prevalence \n(%)")
LA_pastprev_plot
LA_density_plot <- ggplot() +
geom_line(data = LilApp_allpara, aes(x = Julian.Day, y = Density), size = 1.1, alpha = 0.7) +
geom_point(data = LilApp_allpara, aes(x = Julian.Day, y = Density, fill=sampling), size = 5, alpha = 1.0, shape = 21, show.legend = FALSE) +
scale_y_log10() +
scale_x_continuous(breaks = c(213, 244, 274, 305), labels = c("Aug-1", "Sept-1","Oct-1", "Nov-1")) +
scale_fill_manual(values=c("#5445b1", "#749dae", "#f3c483",'#7E7E7E')) +
xlab("Date") +
ylab("*D. dentifera* density \n(per m^2)") +
theme(axis.title.y = ggtext::element_markdown())
LA_density_plot
fielddataplot <- plot_grid(LA_pastprev_plot, LA_density_plot, labels = "auto", ncol = 1, align = "v")
fielddataplot
ggsave(here("figures", "fielddataplot.jpg"), fielddataplot, units = "in", width = 5, height = 7, dpi = 300)
# Load the data for analyzing virulence
LA_infected <- readr::read_csv(here("data/LA_infected.csv"))
## Rows: 267 Columns: 13
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (4): rep, treatment, infection, para_time_point
## dbl (7): clone, block, host_time_point, lifespan, avg, spores_tot, clutches
## date (2): death_date, birth_date
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# Get rid of time 4 (which had only one host clone)
LA_infected <- subset(LA_infected, host_time_point != "4")
# Calculating the proportion that never reproduced
LA_infected <- LA_infected %>%
mutate(anyrepro = if_else(clutches == "0", "0", "1"))
LA_infected$anyrepro <- as.numeric(LA_infected$anyrepro)
# Adding a column for log spores
LA_infected$logspores <- log10(LA_infected$spores_tot)
LA_infected$lnspores <- log(LA_infected$spores_tot)
LA_infected_samplesizes <- LA_infected %>%
count(clone, para_time_point)
# save
write.csv(LA_infected_samplesizes, here("tables", "LA_infected_samplesizes.csv"), row.names = F)
LA_infected$clone <- as.factor(as.character(LA_infected$clone))
LA_infected$host_time_point <- as.factor(as.character(LA_infected$host_time_point))
LA_infected$para_time_point <- as.factor(as.character(LA_infected$para_time_point))
LA_infected$exposed <- ifelse(LA_infected$para_time_point == "none", 'unexposed', 'infected')
# Note: all of the animals in this dataset are either controls or infected, so don't need to worry about exposed but uninfected animals
LA_infected_summary <- LA_infected %>%
group_by(exposed, host_time_point, clone, para_time_point) %>%
summarise(meanlifespan = mean(lifespan), meanclutches = mean(clutches))
## `summarise()` has grouped output by 'exposed', 'host_time_point', 'clone'. You can override using the `.groups` argument.
lifespaninfvcontrol <- ggplot(LA_infected_summary,
aes(x=exposed, y=meanlifespan)) +
geom_violin() +
geom_jitter(shape=16, position=position_jitter(width=0.3, height=0), alpha = 0.5, color = '#5c1a33', size = 2) +
ylab("Mean lifespan \n(days)") +
theme_cowplot() +
theme(axis.title.x = element_blank())
lifespaninfvcontrol
# making a plot comparing exposed vs. unexposed
fecundityinfvcontrol <- ggplot(LA_infected_summary,
aes(x=exposed, y=meanclutches)) +
geom_violin() +
geom_jitter(shape=16, position=position_jitter(width=0.3, height=0), alpha = 0.5, color = '#5c1a33', size = 2) +
ylab("Mean # clutches") +
theme_cowplot() +
theme(axis.title.x = element_blank())
fecundityinfvcontrol
Let’s arrange the lifespan & fecundity plots into a single figure:
infvcontrolplot <- plot_grid(lifespaninfvcontrol, fecundityinfvcontrol, labels = "auto", ncol = 1, align = "v")
infvcontrolplot
ggsave(here("figures", "infvcontrolplot.jpg"), infvcontrolplot, units = "in", width = 5, height = 8, dpi = 300)
infvcontrollifespanmodel<-glmer(lifespan ~ exposed + (1 | clone), family = "poisson", data=LA_infected)
summary(infvcontrollifespanmodel)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: lifespan ~ exposed + (1 | clone)
## Data: LA_infected
##
## AIC BIC logLik deviance df.resid
## 2047.5 2058.2 -1020.8 2041.5 256
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.7830 -0.9877 0.1476 0.9803 4.7174
##
## Random effects:
## Groups Name Variance Std.Dev.
## clone (Intercept) 0.003289 0.05735
## Number of obs: 259, groups: clone, 28
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.53313 0.01679 210.41 <2e-16 ***
## exposedunexposed 0.36282 0.02649 13.69 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## exposdnxpsd -0.353
# Testing for overdispersion (using code from Ben Bolker via Michelle Fearon)
overdisp_fun <- function(model) {
rdf <- df.residual(model)
rp <- residuals(model,type="pearson")
Pearson.chisq <- sum(rp^2)
prat <- Pearson.chisq/rdf
pval <- pchisq(Pearson.chisq, df=rdf, lower.tail=FALSE)
c(chisq=Pearson.chisq,ratio=prat,rdf=rdf,p=pval)
}
overdisp_fun(infvcontrollifespanmodel)
## chisq ratio rdf p
## 5.605749e+02 2.189746e+00 2.560000e+02 7.929284e-25
#very overdispersed, so let's go to a negative binomial distribution
infvcontrollifespanmodel.nb <-glmer.nb(lifespan ~ exposed + (1 | clone), data=LA_infected)
overdisp_fun(infvcontrollifespanmodel.nb)
## chisq ratio rdf p
## 231.3587907 0.9072894 255.0000000 0.8534423
summary(infvcontrollifespanmodel.nb)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(23.7957) ( log )
## Formula: lifespan ~ exposed + (1 | clone)
## Data: LA_infected
##
## AIC BIC logLik deviance df.resid
## 1913.6 1927.9 -952.8 1905.6 255
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.3581 -0.5694 0.0871 0.6945 2.8266
##
## Random effects:
## Groups Name Variance Std.Dev.
## clone (Intercept) 7.208e-05 0.00849
## Number of obs: 259, groups: clone, 28
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.53253 0.01835 192.540 <2e-16 ***
## exposedunexposed 0.36737 0.04151 8.849 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## exposdnxpsd -0.439
infvcontrolrepro<-glmer(clutches ~ exposed + (1 | clone), family = "poisson", data=LA_infected)
overdisp_fun(infvcontrolrepro)
## chisq ratio rdf p
## 8.038865e+02 3.140182e+00 2.560000e+02 7.119246e-58
infvcontrolrepro.nb <-glmer.nb(clutches ~ exposed + (1 | clone), data=LA_infected)
overdisp_fun(infvcontrolrepro.nb)
## chisq ratio rdf p
## 303.84415988 1.19154573 255.00000000 0.01934465
summary(infvcontrolrepro.nb)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(0.8023) ( log )
## Formula: clutches ~ exposed + (1 | clone)
## Data: LA_infected
##
## AIC BIC logLik deviance df.resid
## 836.9 851.1 -414.5 828.9 255
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.8719 -0.5923 -0.4879 0.2711 6.3781
##
## Random effects:
## Groups Name Variance Std.Dev.
## clone (Intercept) 0.466 0.6827
## Number of obs: 259, groups: clone, 28
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.4457 0.1873 -2.379 0.0174 *
## exposedunexposed 2.8101 0.2343 11.995 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## exposdnxpsd -0.281
# The dataset we were working with above only has infected animals or controls. Let's upload the dataset that will also include animals that were exposed but did not become infected
LA_full <- readr::read_csv(here("data/LA_full.csv"))
## Rows: 374 Columns: 13
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr (4): rep, treatment, infection, para_time_point
## dbl (7): clone, block, host_time_point, lifespan, avg, spores_tot, clutches
## date (2): death_date, birth_date
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# since we want to look at infectivity, let's get rid of the control/unexposed animals
LA_exposed <- LA_full %>%
subset(para_time_point != "none")
# get rid of time 4 hosts, since we didn't have parasites from that time point
LA_exposed <- LA_exposed %>%
subset(host_time_point != '4')
# filter out ones that are missing spore counts
LA_exposed <- LA_exposed %>%
drop_na(spores_tot)
# Looking for clones that had very few animals in the experiment
LA_exposed_samplesizes <- LA_exposed %>%
count(clone, para_time_point)
LA_exposed_samplesizes$remove <- ifelse(LA_exposed_samplesizes$n < 3, 'y', 'n')
# save
write.csv(LA_exposed_samplesizes, here("tables", "LA_exposed_samplesizes.csv"), row.names = F)
# Using spore yield to determine if it was infected; this is more sensitive than the 'infection' column, which was based on a quick visual assessment at death -- that inspection missed some that had very low spore densities
LA_exposed$infected <- ifelse(LA_exposed$spores_tot == 0, 0, 1)
# Calculate number infected and total sample sizes
LA_exposed_prevalence <- LA_exposed %>%
group_by(clone, host_time_point, para_time_point) %>%
summarise(total = n(), inf = sum(infected))
## `summarise()` has grouped output by 'clone', 'host_time_point'. You can override using the `.groups` argument.
# Calculate proportion infected per clone
LA_exposed_prevalence$propinf <- LA_exposed_prevalence$inf/LA_exposed_prevalence$total
LA_exposed_prevalence$host_time_point <- as.character(LA_exposed_prevalence$host_time_point)
infprevplot <- ggplot(LA_exposed_prevalence,
aes(x=host_time_point, y=propinf, fill=para_time_point)) +
geom_violin(position=position_dodge(1)) +
geom_jitter(shape=16, position=position_jitter(width=0.3, height=0), alpha = 0.5, show.legend = FALSE) +
scale_fill_manual(values=c("#5445b1", "#749dae", "#f3c483")) +
scale_color_manual(values=c("#5445b1", "#749dae", "#f3c483")) +
xlab("Host time point") +ylab("Proportion infected") +
labs(color= "Parasite time point") + labs(fill="Parasite time point") +
theme_cowplot()
infprevplot
infprevplota <- LA_exposed_prevalence %>%
filter(host_time_point == para_time_point) %>%
ggplot(aes(x=host_time_point, y=propinf, fill=para_time_point)) +
geom_violin(position=position_dodge(1), show.legend = FALSE) +
geom_jitter(shape=16, position=position_jitter(width=0.3, height=0), alpha = 0.5, show.legend = FALSE) +
scale_fill_manual(values=c("#5445b1", "#749dae", "#f3c483")) +
scale_color_manual(values=c("#5445b1", "#749dae", "#f3c483")) +
ylim(0,1.00) +
xlab("Host time point") +ylab("Proportion infected (contemporary parasites)") +
labs(color= "Parasite time point") + labs(fill="Parasite time point") +
theme_cowplot()
infprevplota
infprevplotb <- LA_exposed_prevalence %>%
filter(host_time_point == 3) %>%
ggplot(aes(x=para_time_point, y=propinf, fill=para_time_point)) +
geom_violin(position=position_dodge(1), show.legend = FALSE) +
geom_jitter(shape=16, position=position_jitter(width=0.3, height=0), alpha = 0.5, show.legend = FALSE) +
scale_fill_manual(values=c("#C0C0C0", "#f3c483")) +
scale_color_manual(values=c("#C0C0C0", "#f3c483")) +
ylim(0,1.00) +
xlab("Parasite time point") +ylab("Proportion infected (time 3 hosts)") +
theme_cowplot()
infprevplotb
infprevplot2 <- plot_grid(infprevplota, infprevplotb, labels = "auto", ncol = 2, align = "v")
infprevplot2
ggsave(here("figures", "infprevplot.jpg"), infprevplot2, units = "in", width = 8, height = 5, dpi = 300)
LA_exposed_prevalence_sametime <- subset(LA_exposed_prevalence, host_time_point == para_time_point)
prevmodel1<-glmer(cbind(inf, total-inf) ~ para_time_point + (1 | clone), family='binomial', data=LA_exposed_prevalence_sametime)
# Testing for overdispersion (using code from Ben Bolker via Michelle Fearon)
overdisp_fun <- function(model) {
rdf <- df.residual(model)
rp <- residuals(model,type="pearson")
Pearson.chisq <- sum(rp^2)
prat <- Pearson.chisq/rdf
pval <- pchisq(Pearson.chisq, df=rdf, lower.tail=FALSE)
c(chisq=Pearson.chisq,ratio=prat,rdf=rdf,p=pval)
}
overdisp_fun(prevmodel1)
## chisq ratio rdf p
## 15.3492797 0.6673600 23.0000000 0.8816798
summary(prevmodel1)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(inf, total - inf) ~ para_time_point + (1 | clone)
## Data: LA_exposed_prevalence_sametime
##
## AIC BIC logLik deviance df.resid
## 102.2 107.4 -47.1 94.2 23
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.60835 -0.54861 -0.08292 0.75747 1.13735
##
## Random effects:
## Groups Name Variance Std.Dev.
## clone (Intercept) 0.3203 0.5659
## Number of obs: 27, groups: clone, 27
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.1359 0.4972 0.273 0.785
## para_time_point2 0.4872 0.5938 0.820 0.412
## para_time_point3 0.8869 0.5818 1.524 0.127
##
## Correlation of Fixed Effects:
## (Intr) pr_t_2
## par_tm_pnt2 -0.851
## par_tm_pnt3 -0.873 0.750
LA_exposed_prevalence_sametime$para_time_point <- as.factor(LA_exposed_prevalence_sametime$para_time_point)
# Testing for overall effect of parasite time point
drop1(prevmodel1,test="Chisq")
## Single term deletions
##
## Model:
## cbind(inf, total - inf) ~ para_time_point + (1 | clone)
## npar AIC LRT Pr(Chi)
## <none> 102.17
## para_time_point 2 100.88 2.7111 0.2578
# Now comparing the two sets of spores that were used to infect time 3 hosts
LA_exposed_prevalence_time3 <- subset(LA_exposed_prevalence, host_time_point == '3')
prevmodel2<-glmer(cbind(inf, total-inf) ~ para_time_point + (1 | clone), family='binomial', data=LA_exposed_prevalence_time3)
overdisp_fun(prevmodel2)
## chisq ratio rdf p
## 27.7920278 1.1580012 24.0000000 0.2689096
summary(prevmodel2)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: cbind(inf, total - inf) ~ para_time_point + (1 | clone)
## Data: LA_exposed_prevalence_time3
##
## AIC BIC logLik deviance df.resid
## 94.4 98.3 -44.2 88.4 24
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.3300 -0.8751 0.2081 0.8352 1.4208
##
## Random effects:
## Groups Name Variance Std.Dev.
## clone (Intercept) 0.1248 0.3532
## Number of obs: 27, groups: clone, 14
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.1495 0.2654 4.331 1.49e-05 ***
## para_time_point3 -0.1808 0.3297 -0.548 0.584
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## par_tm_pnt3 -0.673
# Getting chi-square to be parallel to initial analysis
drop1(prevmodel2,test="Chisq")
## Single term deletions
##
## Model:
## cbind(inf, total - inf) ~ para_time_point + (1 | clone)
## npar AIC LRT Pr(Chi)
## <none> 94.402
## para_time_point 1 92.703 0.30031 0.5837
Comparing hosts exposed to parasites from the same time point
LA_sametime <- LA_infected %>%
filter(para_time_point != "none")
LA_sametime$para_time_point <- as.factor(LA_sametime$para_time_point)
LA_sametime$host_time_point <- as.factor(as.character(LA_sametime$host_time_point))
LA_sametime$para_time_point <- droplevels(LA_sametime$para_time_point)
LA_sametime <- subset(LA_sametime, host_time_point == para_time_point)
LA_sametime$paragrowth <- (LA_sametime$spores_tot)/(LA_sametime$lifespan)
LA_sametime_summary <- LA_sametime %>%
group_by(clone, para_time_point) %>%
summarise(meanlifespan = mean(lifespan), meanclutches = mean(clutches), meanspores = mean(spores_tot), meanlogspores = mean(logspores), meanparagrowth = mean(paragrowth))
## `summarise()` has grouped output by 'clone'. You can override using the `.groups` argument.
contemplifespanplot <- ggplot(LA_sametime_summary,
aes(x=para_time_point, y=meanlifespan, fill=para_time_point)) +
geom_violin(show.legend = FALSE) +
geom_jitter(shape=16, position=position_jitter(width=0.3, height=0, seed = 1), alpha = 0.5, show.legend = FALSE) +
scale_fill_manual(values=c("#5445b1", "#749dae", "#f3c483")) +
scale_color_manual(values=c("#5445b1", "#749dae", "#f3c483")) +
xlab("Parasite time point") +ylab("Mean lifespan \n(days)") +
labs(color= "Parasite time point") + labs(fill="Parasite time point") +
theme_cowplot()
contemplifespanplot
contemporarylifespan <-glmer(lifespan ~ host_time_point + (1 | clone), family = "poisson", data=LA_sametime)
overdisp_fun(contemporarylifespan)
## chisq ratio rdf p
## 2.543648e+02 1.829963e+00 1.390000e+02 8.579296e-09
summary(contemporarylifespan)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: lifespan ~ host_time_point + (1 | clone)
## Data: LA_sametime
##
## AIC BIC logLik deviance df.resid
## 1080.7 1092.5 -536.3 1072.7 139
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.0520 -0.9930 -0.0005 0.9847 4.0506
##
## Random effects:
## Groups Name Variance Std.Dev.
## clone (Intercept) 0.009985 0.09992
## Number of obs: 143, groups: clone, 27
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.49051 0.07526 46.376 <2e-16 ***
## host_time_point2 0.09306 0.08537 1.090 0.276
## host_time_point3 0.06512 0.08299 0.785 0.433
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) hst__2
## hst_tm_pnt2 -0.878
## hst_tm_pnt3 -0.910 0.799
contemporarylifespan.nb <-glmer.nb(lifespan ~ host_time_point + (1 | clone), data=LA_sametime)
overdisp_fun(contemporarylifespan.nb)
## chisq ratio rdf p
## 136.4245942 0.9885840 138.0000000 0.5219351
summary(contemporarylifespan.nb)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(28.2922) ( log )
## Formula: lifespan ~ host_time_point + (1 | clone)
## Data: LA_sametime
##
## AIC BIC logLik deviance df.resid
## 1038.2 1053.0 -514.1 1028.2 138
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.9149 -0.7669 0.0226 0.6781 3.6005
##
## Random effects:
## Groups Name Variance Std.Dev.
## clone (Intercept) 0.001694 0.04116
## Number of obs: 143, groups: clone, 27
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.41497 0.07576 45.077 <2e-16 ***
## host_time_point2 0.15898 0.08331 1.908 0.0563 .
## host_time_point3 0.14369 0.08314 1.728 0.0840 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) hst__2
## hst_tm_pnt2 -0.885
## hst_tm_pnt3 -0.924 0.817
drop1(contemporarylifespan.nb,test="Chisq")
## Single term deletions
##
## Model:
## lifespan ~ host_time_point + (1 | clone)
## npar AIC LRT Pr(Chi)
## <none> 1038.2
## host_time_point 2 1037.1 2.9334 0.2307
emmeans(contemporarylifespan.nb, specs = pairwise ~ host_time_point)
## $emmeans
## host_time_point emmean SE df asymp.LCL asymp.UCL
## 1 3.41 0.0758 Inf 3.27 3.56
## 2 3.57 0.0389 Inf 3.50 3.65
## 3 3.56 0.0317 Inf 3.50 3.62
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df z.ratio p.value
## 1 - 2 -0.1590 0.0833 Inf -1.908 0.1363
## 1 - 3 -0.1437 0.0831 Inf -1.728 0.1947
## 2 - 3 0.0153 0.0504 Inf 0.304 0.9504
##
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimates
contempreproplot <- ggplot(LA_sametime_summary,
aes(x=para_time_point, y=meanclutches, fill=para_time_point)) +
geom_violin(show.legend = FALSE) +
geom_jitter(shape=16, position=position_jitter(width=0.3, height=0, seed = 1), alpha = 0.5, show.legend = FALSE) +
scale_fill_manual(values=c("#5445b1", "#749dae", "#f3c483")) +
scale_color_manual(values=c("#5445b1", "#749dae", "#f3c483")) +
xlab("Parasite time point") +ylab("Mean number of clutches") +
labs(color= "Parasite time point") + labs(fill="Parasite time point") +
theme_cowplot()
contempreproplot
contemporaryclutches <-glmer(clutches ~ host_time_point + (1 | clone), family = "poisson", data=LA_sametime)
overdisp_fun(contemporaryclutches)
## chisq ratio rdf p
## 2.724757e+02 1.960257e+00 1.390000e+02 1.038315e-10
summary(contemporaryclutches)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: clutches ~ host_time_point + (1 | clone)
## Data: LA_sametime
##
## AIC BIC logLik deviance df.resid
## 445.8 457.6 -218.9 437.8 139
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.2287 -0.7444 -0.4875 0.0935 4.9567
##
## Random effects:
## Groups Name Variance Std.Dev.
## clone (Intercept) 1.627 1.275
## Number of obs: 143, groups: clone, 27
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.1837 0.7839 -0.234 0.815
## host_time_point2 -0.1733 0.8959 -0.193 0.847
## host_time_point3 -0.7804 0.8713 -0.896 0.370
##
## Correlation of Fixed Effects:
## (Intr) hst__2
## hst_tm_pnt2 -0.858
## hst_tm_pnt3 -0.875 0.760
contemporaryclutches.nb <-glmer.nb(clutches ~ host_time_point + (1 | clone), data=LA_sametime)
overdisp_fun(contemporaryclutches.nb)
## chisq ratio rdf p
## 106.5188352 0.7718756 138.0000000 0.9783568
summary(contemporaryclutches.nb)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(0.3527) ( log )
## Formula: clutches ~ host_time_point + (1 | clone)
## Data: LA_sametime
##
## AIC BIC logLik deviance df.resid
## 344.8 359.6 -167.4 334.8 138
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.5582 -0.4579 -0.4028 0.1347 3.7312
##
## Random effects:
## Groups Name Variance Std.Dev.
## clone (Intercept) 0.8264 0.9091
## Number of obs: 143, groups: clone, 27
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.2162 0.7690 0.281 0.779
## host_time_point2 -0.4843 0.8720 -0.555 0.579
## host_time_point3 -1.0601 0.8472 -1.251 0.211
##
## Correlation of Fixed Effects:
## (Intr) hst__2
## hst_tm_pnt2 -0.863
## hst_tm_pnt3 -0.884 0.784
drop1(contemporaryclutches.nb,test="Chisq")
## Single term deletions
##
## Model:
## clutches ~ host_time_point + (1 | clone)
## npar AIC LRT Pr(Chi)
## <none> 344.82
## host_time_point 2 342.78 1.9565 0.376
emmeans(contemporaryclutches.nb, specs = pairwise ~ host_time_point)
## $emmeans
## host_time_point emmean SE df asymp.LCL asymp.UCL
## 1 0.216 0.769 Inf -1.29 1.7234
## 2 -0.268 0.441 Inf -1.13 0.5965
## 3 -0.844 0.397 Inf -1.62 -0.0658
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df z.ratio p.value
## 1 - 2 0.484 0.872 Inf 0.555 0.8437
## 1 - 3 1.060 0.847 Inf 1.251 0.4229
## 2 - 3 0.576 0.566 Inf 1.017 0.5658
##
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimates
contempsporesplot <- ggplot(LA_sametime_summary,
aes(x=para_time_point, y=meanspores, fill=para_time_point)) +
geom_violin(show.legend = FALSE) +
geom_jitter(shape=16, position=position_jitter(width=0.3, height=0, seed = 1), alpha = 0.5, show.legend = FALSE) +
scale_fill_manual(values=c("#5445b1", "#749dae", "#f3c483")) +
scale_color_manual(values=c("#5445b1", "#749dae", "#f3c483")) +
xlab("Parasite time point") +ylab("Mean number of \nspores per infected host") +
labs(color= "Parasite time point") + labs(fill="Parasite time point") +
scale_y_log10() +
theme_cowplot()
contempsporesplot
contemporarylogspores <- lmer(logspores ~ host_time_point + (1 | clone), data=LA_sametime)
overdisp_fun(contemporarylogspores)
## chisq ratio rdf p
## 48.7586242 0.3533234 138.0000000 1.0000000
summary(contemporarylogspores)
## Linear mixed model fit by REML ['lmerMod']
## Formula: logspores ~ host_time_point + (1 | clone)
## Data: LA_sametime
##
## REML criterion at convergence: 284.6
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.8804 -0.3749 0.1072 0.6637 1.5019
##
## Random effects:
## Groups Name Variance Std.Dev.
## clone (Intercept) 0.05924 0.2434
## Residual 0.37508 0.6124
## Number of obs: 143, groups: clone, 27
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 4.96189 0.21368 23.221
## host_time_point2 -0.07965 0.24515 -0.325
## host_time_point3 -0.25471 0.23558 -1.081
##
## Correlation of Fixed Effects:
## (Intr) hst__2
## hst_tm_pnt2 -0.872
## hst_tm_pnt3 -0.907 0.791
drop1(contemporarylogspores,test="Chisq")
## Single term deletions
##
## Model:
## logspores ~ host_time_point + (1 | clone)
## npar AIC LRT Pr(Chi)
## <none> 287.95
## host_time_point 2 286.11 2.166 0.3386
emmeans(contemporarylogspores, specs = pairwise ~ host_time_point)
## $emmeans
## host_time_point emmean SE df lower.CL upper.CL
## 1 4.96 0.2201 18.9 4.50 5.42
## 2 4.88 0.1209 22.2 4.63 5.13
## 3 4.71 0.0995 20.1 4.50 4.91
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## 1 - 2 0.0796 0.251 19.6 0.317 0.9462
## 1 - 3 0.2547 0.242 19.1 1.055 0.5527
## 2 - 3 0.1751 0.157 21.3 1.118 0.5136
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimates
#R3 asked about the non-log-transformed data, so let's do that
contemporaryspores <- lmer(spores_tot ~ host_time_point + (1 | clone), data=LA_sametime)
overdisp_fun(contemporaryspores)
## chisq ratio rdf p
## 3.097726e+12 2.244729e+10 1.380000e+02 0.000000e+00
summary(contemporaryspores)
## Linear mixed model fit by REML ['lmerMod']
## Formula: spores_tot ~ host_time_point + (1 | clone)
## Data: LA_sametime
##
## REML criterion at convergence: 3760.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.4286 -0.5585 -0.3070 0.3365 4.7327
##
## Random effects:
## Groups Name Variance Std.Dev.
## clone (Intercept) 2.312e+09 48086
## Residual 2.339e+10 152932
## Number of obs: 143, groups: clone, 27
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 136891 47944 2.855
## host_time_point2 35071 55183 0.636
## host_time_point3 -22837 52941 -0.431
##
## Correlation of Fixed Effects:
## (Intr) hst__2
## hst_tm_pnt2 -0.869
## hst_tm_pnt3 -0.906 0.787
drop1(contemporaryspores,test="Chisq")
## Single term deletions
##
## Model:
## spores_tot ~ host_time_point + (1 | clone)
## npar AIC LRT Pr(Chi)
## <none> 3837.7
## host_time_point 2 3836.6 2.914 0.2329
contempparagrowthplot <- ggplot(LA_sametime_summary,
aes(x=para_time_point, y=meanparagrowth, fill=para_time_point)) +
geom_violin(show.legend = FALSE) +
geom_jitter(shape=16, position=position_jitter(width=0.3, height=0, seed = 1), alpha = 0.5, show.legend = FALSE) +
scale_fill_manual(values=c("#5445b1", "#749dae", "#f3c483")) +
scale_color_manual(values=c("#5445b1", "#749dae", "#f3c483")) +
xlab("Parasite time point") +ylab("Parasite growth rate \n(spores per day)") +
labs(color= "Parasite time point") + labs(fill="Parasite time point") +
theme_cowplot()
contempparagrowthplot
contemporaryparasitegrowth <- lmer(paragrowth ~ host_time_point + (1 | clone), data=LA_sametime)
overdisp_fun(contemporaryparasitegrowth)
## chisq ratio rdf p
## 1994481979 14452768 138 0
summary(contemporaryparasitegrowth)
## Linear mixed model fit by REML ['lmerMod']
## Formula: paragrowth ~ host_time_point + (1 | clone)
## Data: LA_sametime
##
## REML criterion at convergence: 2736
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6309 -0.5728 -0.2677 0.3091 4.2680
##
## Random effects:
## Groups Name Variance Std.Dev.
## clone (Intercept) 2065795 1437
## Residual 15245387 3905
## Number of obs: 143, groups: clone, 27
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 4279.9 1313.2 3.259
## host_time_point2 538.4 1507.8 0.357
## host_time_point3 -1219.3 1448.3 -0.842
##
## Correlation of Fixed Effects:
## (Intr) hst__2
## hst_tm_pnt2 -0.871
## hst_tm_pnt3 -0.907 0.790
drop1(contemporaryparasitegrowth,test="Chisq")
## Single term deletions
##
## Model:
## paragrowth ~ host_time_point + (1 | clone)
## npar AIC LRT Pr(Chi)
## <none> 2791.7
## host_time_point 2 2791.4 3.6498 0.1612
emmeans(contemporaryparasitegrowth, specs = pairwise ~ host_time_point)
## $emmeans
## host_time_point emmean SE df lower.CL upper.CL
## 1 4280 1357 17.4 1422 7138
## 2 4818 746 22.0 3272 6365
## 3 3061 613 19.9 1782 4339
##
## Degrees-of-freedom method: kenward-roger
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df t.ratio p.value
## 1 - 2 -538 1548 18.4 -0.348 0.9358
## 1 - 3 1219 1489 17.8 0.819 0.6965
## 2 - 3 1758 965 21.1 1.821 0.1870
##
## Degrees-of-freedom method: kenward-roger
## P value adjustment: tukey method for comparing a family of 3 estimates
Comparing hosts from time 3 that were exposed to parasites from time 1 vs. time 3
LA_paraevol <- LA_infected %>%
filter(host_time_point == "3") %>%
filter(para_time_point != "none")
LA_paraevol$paragrowth <- (LA_paraevol$spores_tot)/(LA_paraevol$lifespan)
LA_paraevol_summary <- LA_paraevol %>%
group_by(clone, para_time_point) %>%
summarise(meanlifespan = mean(lifespan), meanclutches = mean(clutches), meanspores = mean(spores_tot), meanlogspores = mean(logspores), meanparagrowth = mean(paragrowth))
## `summarise()` has grouped output by 'clone'. You can override using the `.groups` argument.
lifespanparaevol <- ggplot(LA_paraevol_summary,
aes(x=para_time_point, y=meanlifespan, fill=para_time_point)) +
geom_violin(show.legend = FALSE) +
geom_jitter(shape=16, position=position_jitter(width=0.3, height=0, seed = 1), alpha = 0.5, show.legend = FALSE) +
scale_fill_manual(values=c("#C0C0C0", "#f3c483")) +
scale_color_manual(values=c("#C0C0C0", "#f3c483")) +
xlab("Parasite time point") +ylab("Mean lifespan \n(days)") +
labs(color= "Parasite time point") + labs(fill="Parasite time point") +
theme_cowplot()
lifespanparaevol
paraevollifespan <-glmer(lifespan ~ para_time_point + (1 | clone), family = "poisson", data=LA_paraevol)
overdisp_fun(paraevollifespan)
## chisq ratio rdf p
## 2.459747e+02 1.732216e+00 1.420000e+02 1.428517e-07
summary(paraevollifespan)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: lifespan ~ para_time_point + (1 | clone)
## Data: LA_paraevol
##
## AIC BIC logLik deviance df.resid
## 1057.3 1066.2 -525.7 1051.3 142
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.2121 -0.8940 0.1315 0.8178 3.1930
##
## Random effects:
## Groups Name Variance Std.Dev.
## clone (Intercept) 0.003669 0.06057
## Number of obs: 145, groups: clone, 14
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.50151 0.02748 127.421 <2e-16 ***
## para_time_point3 0.05260 0.02914 1.805 0.0711 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## par_tm_pnt3 -0.577
paraevollifespan.nb <-glmer.nb(lifespan ~ para_time_point + (1 | clone), data=LA_paraevol)
overdisp_fun(paraevollifespan.nb)
## chisq ratio rdf p
## 137.3167343 0.9738775 141.0000000 0.5720026
summary(paraevollifespan.nb)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(38.342) ( log )
## Formula: lifespan ~ para_time_point + (1 | clone)
## Data: LA_paraevol
##
## AIC BIC logLik deviance df.resid
## 1027.4 1039.3 -509.7 1019.4 141
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -3.1564 -0.6213 0.0274 0.6221 2.4806
##
## Random effects:
## Groups Name Variance Std.Dev.
## clone (Intercept) 0.001195 0.03458
## Number of obs: 145, groups: clone, 14
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.50845 0.03044 115.264 <2e-16 ***
## para_time_point3 0.04913 0.03937 1.248 0.212
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## par_tm_pnt3 -0.680
reproparaevol <- ggplot(LA_paraevol_summary,
aes(x=para_time_point, y=meanclutches, fill=para_time_point)) +
geom_violin(show.legend = FALSE) +
geom_jitter(shape=16, position=position_jitter(width=0.3, height=0, seed = 1), alpha = 0.5, show.legend = FALSE) +
scale_fill_manual(values=c("#C0C0C0", "#f3c483")) +
scale_color_manual(values=c("#C0C0C0", "#f3c483")) +
xlab("Parasite time point") +ylab("Mean number of clutches") +
labs(color= "Parasite time point") + labs(fill="Parasite time point") +
theme_cowplot()
reproparaevol
paraevolclutches <-glmer(clutches ~ para_time_point + (1 | clone), family = "poisson", data=LA_paraevol)
overdisp_fun(paraevolclutches)
## chisq ratio rdf p
## 3.081137e+02 2.169815e+00 1.420000e+02 2.583562e-14
summary(paraevolclutches)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: poisson ( log )
## Formula: clutches ~ para_time_point + (1 | clone)
## Data: LA_paraevol
##
## AIC BIC logLik deviance df.resid
## 344.7 353.6 -169.3 338.7 142
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.5607 -0.7308 -0.5406 -0.2792 7.3191
##
## Random effects:
## Groups Name Variance Std.Dev.
## clone (Intercept) 0.5547 0.7448
## Number of obs: 145, groups: clone, 14
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.0469 0.2896 -3.615 0.000301 ***
## para_time_point3 0.4432 0.2380 1.862 0.062623 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## par_tm_pnt3 -0.516
paraevolclutches.nb <-glmer.nb(clutches ~ para_time_point + (1 | clone), data=LA_paraevol)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00817007 (tol = 0.002, component 1)
## boundary (singular) fit: see ?isSingular
overdisp_fun(paraevolclutches.nb)
## chisq ratio rdf p
## 158.6689373 1.1253116 141.0000000 0.1467461
summary(paraevolclutches.nb)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(0.2552) ( log )
## Formula: clutches ~ para_time_point + (1 | clone)
## Data: LA_paraevol
##
## AIC BIC logLik deviance df.resid
## 276.5 288.4 -134.3 268.5 141
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.4292 -0.4292 -0.3989 0.2190 7.3489
##
## Random effects:
## Groups Name Variance Std.Dev.
## clone (Intercept) 2.695e-10 1.642e-05
## Number of obs: 145, groups: clone, 14
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.8615 0.2975 -2.895 0.00379 **
## para_time_point3 0.4492 0.4023 1.117 0.26419
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr)
## par_tm_pnt3 -0.739
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
sporesparaevol <- ggplot(LA_paraevol_summary,
aes(x=para_time_point, y=meanspores, fill=para_time_point)) +
geom_violin(show.legend = FALSE) +
geom_jitter(shape=16, position=position_jitter(width=0.3, height=0, seed = 1), alpha = 0.5, show.legend = FALSE) +
scale_fill_manual(values=c("#C0C0C0", "#f3c483")) +
scale_color_manual(values=c("#C0C0C0", "#f3c483")) +
xlab("Parasite time point") +ylab("Mean number of \nspores per infected host") +
labs(color= "Parasite time point") + labs(fill="Parasite time point") +
scale_y_log10() +
theme_cowplot()
sporesparaevol
paraevollogspores <- lmer(logspores ~ para_time_point + (1 | clone), data=LA_paraevol)
overdisp_fun(paraevollogspores)
## chisq ratio rdf p
## 44.3092364 0.3142499 141.0000000 1.0000000
summary(paraevollogspores)
## Linear mixed model fit by REML ['lmerMod']
## Formula: logspores ~ para_time_point + (1 | clone)
## Data: LA_paraevol
##
## REML criterion at convergence: 259.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -4.1494 -0.3281 0.1027 0.6344 1.7752
##
## Random effects:
## Groups Name Variance Std.Dev.
## clone (Intercept) 0.02438 0.1562
## Residual 0.32184 0.5673
## Number of obs: 145, groups: clone, 14
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 5.02231 0.08202 61.233
## para_time_point3 -0.31407 0.09616 -3.266
##
## Correlation of Fixed Effects:
## (Intr)
## par_tm_pnt3 -0.613
drop1(paraevollogspores,test="Chisq")
## Single term deletions
##
## Model:
## logspores ~ para_time_point + (1 | clone)
## npar AIC LRT Pr(Chi)
## <none> 260.88
## para_time_point 1 269.32 10.444 0.00123 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#R3 asked about the analysis if we don't log transform the data
paraevolspores <- lmer(spores_tot ~ para_time_point + (1 | clone), data=LA_paraevol)
overdisp_fun(paraevolspores)
## chisq ratio rdf p
## 2.851456e+12 2.022309e+10 1.410000e+02 0.000000e+00
summary(paraevolspores)
## Linear mixed model fit by REML ['lmerMod']
## Formula: spores_tot ~ para_time_point + (1 | clone)
## Data: LA_paraevol
##
## REML criterion at convergence: 3813.2
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.1691 -0.5994 -0.3277 0.3276 4.7299
##
## Random effects:
## Groups Name Variance Std.Dev.
## clone (Intercept) 7.469e+08 27330
## Residual 2.042e+10 142903
## Number of obs: 145, groups: clone, 14
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 168349 18911 8.902
## para_time_point3 -54543 24030 -2.270
##
## Correlation of Fixed Effects:
## (Intr)
## par_tm_pnt3 -0.658
drop1(paraevolspores,test="Chisq")
## Single term deletions
##
## Model:
## spores_tot ~ para_time_point + (1 | clone)
## npar AIC LRT Pr(Chi)
## <none> 3864.1
## para_time_point 1 3867.3 5.157 0.02315 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
paragrowthparaevol <- ggplot(LA_paraevol_summary,
aes(x=para_time_point, y=meanparagrowth, fill=para_time_point)) +
geom_violin(show.legend = FALSE) +
geom_jitter(shape=16, position=position_jitter(width=0.3, height=0, seed = 1), alpha = 0.5, show.legend = FALSE) +
scale_fill_manual(values=c("#C0C0C0", "#f3c483")) +
scale_color_manual(values=c("#C0C0C0", "#f3c483")) +
xlab("Parasite time point") +ylab("Parasite growth rate \n(spores per day)") +
labs(color= "Parasite time point") + labs(fill="Parasite time point") +
theme_cowplot()
paragrowthparaevol
paraevolparasitegrowth <- lmer(paragrowth ~ para_time_point + (1 | clone), data=LA_paraevol)
summary(paraevolparasitegrowth)
## Linear mixed model fit by REML ['lmerMod']
## Formula: paragrowth ~ para_time_point + (1 | clone)
## Data: LA_paraevol
##
## REML criterion at convergence: 2778.4
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.3042 -0.6355 -0.2704 0.3463 4.5048
##
## Random effects:
## Groups Name Variance Std.Dev.
## clone (Intercept) 697521 835.2
## Residual 14604033 3821.5
## Number of obs: 145, groups: clone, 14
##
## Fixed effects:
## Estimate Std. Error t value
## (Intercept) 5003.8 519.8 9.626
## para_time_point3 -1942.4 644.3 -3.015
##
## Correlation of Fixed Effects:
## (Intr)
## par_tm_pnt3 -0.644
drop1(paraevolparasitegrowth,test="Chisq")
## Single term deletions
##
## Model:
## paragrowth ~ para_time_point + (1 | clone)
## npar AIC LRT Pr(Chi)
## <none> 2815.0
## para_time_point 1 2821.9 8.9619 0.002757 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
eightpanelplot <- plot_grid(contemplifespanplot, lifespanparaevol, contempreproplot, reproparaevol, contempsporesplot, sporesparaevol, contempparagrowthplot, paragrowthparaevol, labels = "auto", ncol = 2, align = "v")
eightpanelplot
# now add the title
title1 <- ggplot() + # Draw ggplot2 plot with text only
annotate("text",
x = 0,
y = 1,
size = 6,
label = "Contemporary hosts",
# hjust = 0.64,
hjust = 0.5,
fontface = "bold") +
theme_void()
title2 <- ggplot() + # Draw ggplot2 plot with text only
annotate("text",
x = 0,
y = 1,
size = 6,
label = "Time 3 hosts",
hjust = 0.5,
fontface = "bold") +
theme_void()
titles <- plot_grid(title1, title2, nrow = 1, hjust = 0)
titles
eightpanelplot_withtitles <- plot_grid(
titles, eightpanelplot,
ncol = 1,
# rel_heights values control vertical title margins
rel_heights = c(0.05, 1)
)
eightpanelplot_withtitles
ggsave(here("figures", "eightpanelplot_withtitles.jpg"), eightpanelplot_withtitles, units = "in", width = 8, height = 12, dpi = 300)
Reviewer 3 was asking about exposed but uninfected animals. Let’s make a new dataframe to look at that.
# get rid of time 4 hosts from the full dataset, since we didn't have parasites from that time point
LA_full <- LA_full %>%
subset(host_time_point != '4')
# filter out ones that are missing spore counts
LA_full <- LA_full %>%
drop_na(spores_tot)
# Looking for clones that had very few animals in the experiment
LA_full_samplesizes <- LA_full %>%
count(clone, para_time_point)
# save
write.csv(LA_full_samplesizes, here("tables", "LA_full_samplesizes.csv"), row.names = F)
# Using spore yield to determine if it was infected; this is more sensitive than the 'infection' column, which was based on a quick visual assessment at death -- that inspection missed some that had very low spore densities
LA_full$infected <- ifelse(LA_full$spores_tot == 0, 0, 1)
LA_full$exposureclass <- ifelse(LA_full$para_time_point == "none", 'control',
ifelse(LA_full$spores_tot > 0, 'infected', 'exposed'))
LA_full_summary <- LA_full %>%
group_by(exposureclass, host_time_point, clone, para_time_point) %>%
summarise(meanlifespan = mean(lifespan), meanclutches = mean(clutches))
## `summarise()` has grouped output by 'exposureclass', 'host_time_point', 'clone'. You can override using the `.groups` argument.
lifespan_exposureclass <- ggplot(LA_full_summary,
aes(x=exposureclass, y=meanlifespan)) +
geom_violin() +
geom_jitter(shape=16, position=position_jitter(width=0.3, height=0), alpha = 0.5, color = '#5c1a33', size = 2) +
ylab("Mean lifespan \n(days)") +
theme_cowplot() +
theme(axis.title.x = element_blank())
lifespan_exposureclass
# making a plot comparing exposed vs. unexposed
repro_exposureclass <- ggplot(LA_full_summary,
aes(x=exposureclass, y=meanclutches)) +
geom_violin() +
geom_jitter(shape=16, position=position_jitter(width=0.3, height=0), alpha = 0.5, color = '#5c1a33', size = 2) +
ylab("Mean number of clutches") +
theme_cowplot() +
theme(axis.title.x = element_blank())
repro_exposureclass
Let’s arrange the lifespan & fecundity plots into a single figure:
exposureclassplot <- plot_grid(lifespan_exposureclass, repro_exposureclass, labels = "auto", ncol = 1, align = "v")
exposureclassplot
ggsave(here("figures", "exposureclassplot.jpg"), exposureclassplot, units = "in", width = 5, height = 8, dpi = 300)
lifespanmodel_3classes.nb <-glmer.nb(lifespan ~ exposureclass + (1 | clone), data=LA_full)
overdisp_fun(lifespanmodel_3classes.nb)
## chisq ratio rdf p
## 269.8401318 0.7754027 348.0000000 0.9992965
summary(lifespanmodel_3classes.nb)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(12.3054) ( log )
## Formula: lifespan ~ exposureclass + (1 | clone)
## Data: LA_full
##
## AIC BIC logLik deviance df.resid
## 2796.9 2816.3 -1393.5 2786.9 348
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -2.6273 -0.4619 0.1239 0.6483 2.3959
##
## Random effects:
## Groups Name Variance Std.Dev.
## clone (Intercept) 0.001463 0.03825
## Number of obs: 353, groups: clone, 28
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.89495 0.04953 78.644 < 2e-16 ***
## exposureclassexposed -0.10544 0.05889 -1.790 0.0734 .
## exposureclassinfected -0.36243 0.05412 -6.696 2.14e-11 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) expsrclssx
## expsrclssxp -0.811
## expsrclssnf -0.894 0.742
emmeans(lifespanmodel_3classes.nb, specs = pairwise ~ exposureclass)
## $emmeans
## exposureclass emmean SE df asymp.LCL asymp.UCL
## control 3.89 0.0495 Inf 3.80 3.99
## exposed 3.79 0.0345 Inf 3.72 3.86
## infected 3.53 0.0242 Inf 3.49 3.58
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df z.ratio p.value
## control - exposed 0.105 0.0589 Inf 1.790 0.1727
## control - infected 0.362 0.0541 Inf 6.696 <.0001
## exposed - infected 0.257 0.0408 Inf 6.292 <.0001
##
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimates
repromodel_3classes.nb <-glmer.nb(clutches ~ exposureclass + (1 | clone), data=LA_full)
summary(repromodel_3classes.nb)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Negative Binomial(1.0935) ( log )
## Formula: clutches ~ exposureclass + (1 | clone)
## Data: LA_full
##
## AIC BIC logLik deviance df.resid
## 1492.3 1511.7 -741.2 1482.3 348
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.0094 -0.6569 -0.5309 0.2700 10.3234
##
## Random effects:
## Groups Name Variance Std.Dev.
## clone (Intercept) 0.1559 0.3948
## Number of obs: 353, groups: clone, 28
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.4587 0.1879 13.083 <2e-16 ***
## exposureclassexposed -0.1459 0.1994 -0.732 0.464
## exposureclassinfected -2.7900 0.1939 -14.391 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) expsrclssx
## expsrclssxp -0.760
## expsrclssnf -0.756 0.698
emmeans(repromodel_3classes.nb, specs = pairwise ~ exposureclass)
## $emmeans
## exposureclass emmean SE df asymp.LCL asymp.UCL
## control 2.459 0.188 Inf 2.090 2.8270
## exposed 2.313 0.135 Inf 2.049 2.5767
## infected -0.331 0.133 Inf -0.593 -0.0698
##
## Results are given on the log (not the response) scale.
## Confidence level used: 0.95
##
## $contrasts
## contrast estimate SE df z.ratio p.value
## control - exposed 0.146 0.199 Inf 0.732 0.7445
## control - infected 2.790 0.194 Inf 14.391 <.0001
## exposed - infected 2.644 0.153 Inf 17.283 <.0001
##
## Results are given on the log (not the response) scale.
## P value adjustment: tukey method for comparing a family of 3 estimates